function I = draw_speckle_pattern(row, col, centers, R)

X = centers(:,1);  % X coordinates of speckle centers
Y = centers(:,2);  % Y coordinates of speckle centers
num = numel(X);    % Number of speckles
A = 255;  % Amplitude of the Gaussian distribution

% Generate coordinate grids for the entire image
[x_grid, y_grid] = meshgrid(1:col, 1:row);

% Calculate the difference in x and y coordinates between each pixel and each speckle center
x_diff = x_grid - reshape(X, 1, 1, num);  % Difference in x-coordinates
y_diff = y_grid - reshape(Y, 1, 1, num);  % Difference in y-coordinates

% Calculate the Gaussian values for all pixels and speckles at once
gauss_values = A * exp(-(x_diff.^2 + y_diff.^2) / (R^2));

% Sum up the Gaussian contributions from all speckles for each pixel
I = sum(gauss_values, 3);  % Sum along the rows to get the total intensity for each pixel

% Ensure pixel values are within the valid range [0, 255]
I(I > 255) = 255;
I = uint8(I);  % Convert to uint8 for proper image representation
end